% clear;
% load('datastatic')
for j=1:length(ss)-1
    s=ss(j);
betas=beta0+beta1*s;
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);
ix=ixx(j);
x=xx(j);
errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j);   
end
errorode(j+1)=0;


for j=1:length(ss)-1
%bs(j)=bs(j+1)-dbs(j+1)*ds;
s=ss(j);
betas=beta0+beta1*s;
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);
ix=ixx(j);
x=xx(j);

for kkj=1:20
db=(bs(j+1)-bs(j))/ds;



errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j);
while abs(errorode(j))>10^(-12)
derrorode=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi-1)*(1-1/psi)*(A-ix-x)/bs(j)^2*(-1))+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*((-1)/ds/bs(j)+(bs(j+1)-bs(j))/ds/bs(j)^2*(-1));
bs(j)=bs(j)-errorode(j)/derrorode;
errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j);   
end





end
ixx(j)=ix;
xx(j)=x;


end








for j=1:length(ss)-1
dbsnu(j)=(bs(j+1)-bs(j))/ds;
end
dbsnu(j+1)=2*dbsnu(j)-dbsnu(j-1);

for j=1:length(ss)
dbs(j)=(A-ix-x)^(-1/psi)*rho/(1-etaX*x/s)*bx(j)^(1/psi);
end

i=1;
while dbsnu(i)>dbs(i) && i<length(ss)
i=i+1;
end
flagstar=i;
ss(flagstar)


% 
% 
% 
% 
% subplot(3,2,1)
% plot(ss,bs,'LineWidth',1.5)
% hold on;
% 
% plot(ss,bs0,'k:','LineWidth',1.5)
% hold on;
% 
% plot(ss,bx,'r--','LineWidth',1.5)
% hold on;
% 
% 
% 
% 
% subplot(3,2,2)
% plot(ss,ixx,'LineWidth',1.5)
% hold on;
% 
% plot(ss,ixx0,'k:','LineWidth',1.5)
% hold on;
% 
% subplot(3,2,3)
% plot(ss,xx,'LineWidth',1.5)
% hold on;
% 
% plot(ss,xx0,'k:','LineWidth',1.5)
% hold on;
% 
% 
% subplot(3,2,4)
% plot(ss,errorode,'LineWidth',1.5)
% hold on;
% 
% 
% subplot(3,2,5)
% plot(ss,ixx-ixx0,'LineWidth',1.5)
% hold on;
% 
% subplot(3,2,6)
% plot(ss,bs-bs0,'LineWidth',1.5)
% hold on;
% 
% % subplot(3,2,5)
% % plot(ss,dbs-dbsan,'LineWidth',1.5)
% % hold on;
